Chapter Three 


State Space Approach 


The state space approach has been introduced in Section 1.3. Due to its 
fundamental importance for control systems, the state space technique will be 
considered thoroughly in this chapter. Both continuous- and discrete-time linear 
time invariant systems will be presented. It has already been pointed out that 
the state space technique represents the modern approach to control system 
theory and its applications. The state space approach is very convenient for 
representation of high-order dimensional and complex systems, and extremely 
efficient for numerical calculations since many efficient and reliable numerical 
algorithms developed in mathematics, especially within the area of numerical 
lineal - algebra, can be used directly. In addition, the state space form is the 
basis for introducing system controllability and observability concepts and many 
modern control theory techniques. 

The state space model of a continuous-time linear system is represented by 
a system of n linear differential equations. In matrix form, it is given by 

^x(t)=x(t) = Ax(t) + Bu(t), x(0) = x o (3.1) 

y (t) = Cx(i) + Du(l) (3.2) 

where x 6 u 6 and y 6 3? p are, respectively, vectors of system 
states, control inputs, and system outputs. The matrix A nXn describes the 
internal behavior of the system, while matrices B" x, \ C pX ", and D pXr represent 
connections between the external world and the system. If there are no direct 
paths between inputs and outputs, which is often the case, the matrix D pXm is 
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zero. It is assumed in this book that all matrices in (3.1) and (3.2) are time 
invariant. Studying linear control systems with time varying coefficient matrices 
requires knowledge of some advanced topics in mathematics (see for example 
Chen, 1984; see also Section 10.1). 

The state space model for linear discrete-time control systems has exactly the 
same form as (3.1) and (3.2) with differential equations replaced by difference 
equations, that is 

x(/r + 1) = A d x.(k) + B d u(k), x(0) = x 0 (3.3) 

y(k) = C d x.(k) + D d u(k) (3.4) 

All vectors and matrices defined in (3.3) and (3.4) have the same dimensions as 
corresponding ones given in (3.1) and (3.2). In this chapter, we present and derive 
in detail the main state space concepts for continuous-time linear control systems 
and then give the corresponding interpretations in the discrete-time domain. 

The chapter is organized as follows. In Section 3.1 several systematic 
methods for obtaining the state space form from differential equations and transfer 
functions are developed. The time response of linear systems given in the state 
form is considered in Section 3.2. The corresponding results for discrete-time 
systems, and the procedure for discretization of continuous-time systems leading 
to discrete-time models, are given in Section 3.3. The concepts of the system 
characteristic equation, eigenvalues, and eigenvectors and their use in control 
system theory are presented in Section 3.4. At the end of the chapter, in Section 
3.5, three MATLAB laboratory experiments are outlined. 

Chapter Objectives 

The dynamical systems considered in this book are either described by 
differential/difference equations or given in the form of system transfer functions. 
One of the goals is to present procedures for obtaining the state space forms either 
from differential/difference equations or from transfer functions. In that respect 
students will be exposed to four standard state space forms, known as canonical 
forms: the phase variable form or controller form, the observer form, the modal 
form, and the Jordan form. 

Another important objective is to show students how to analyze linear 
systems given in the state space form, i.e. how to find responses (state variables 
and outputs) of the corresponding state space models to any input signal (control 
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input). A working knowledge of undergraduate linear algebra and the basic 
theory of differential equations is helpful for complete understanding of this 
chapter. Some useful results on linear algebra arc given in Appendix C. Students 
without a strong background in these topics may consult any undergraduate text, 
(for example, Fraleigh and Beauregard, 1990; Boyce and DiPrima, 1992). 


3.1 State Space Models 

In this section we study state space models of continuous-time linear systems. 
The corresponding results for discrete-time systems, obtained via duality with the 
continuous-time models, are given in Section 3.3. 

The state space model of a continuous-time dynamic system can be derived 
either from the system model given in the time domain by a differential equation 
or from its transfer function representation. Both cases will be considered in this 
section. Four state space forms—the phase variable form (controller form), the 
observer form, the modal form, and the Jordan form—which are often used in 
modern control theory and practice, are presented. 


3.1.1 The State Space Model and Differential Equations 

Consider a general //th-order model of a dynamic system represented by an nth- 
order differential equation 


d n y(t) , 

dt n " r 


& 71 — 1 


d n 1 y(t) 
dt 71 - 1 


+ ••• + 01 


dy(t) 

dt 


+ a 0 y(t) 


= b 


n 


d n u(t) 

dt n 


+ b n -i 


dt 71 - 1 



+ b 0 u(t) 


(3.5) 


At this point we assume that all initial conditions for the above differential 
equation, i.e. y(0~),dy(0~)/dt,...,d n ~ 1 y(0~)/dt n ~ 1 , are equal to zero. We 
will show later how to take into account the effect of initial conditions. 

In order to derive a systematic procedure that transforms a differential 
equation of order n to a state space form representing a system of n first-order 
differential equations, we first start with a simplified version of (3.5), namely we 
study the case when no derivatives with respect to the input are present 


d n y(t) 
dt n ^ 


ttn—l 


d^yjt) 

dt 17 - 1 


CL\ 


dy(t) 


aoy(t) 


u(t) 


(3.6) 


dt 
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Introduce the following (easy to remember) change of variables 

xi (t) = y(t) 
dy(t) 


X-2 (t) = 

X 3 (t) = 


dt 

d 2 y(t) 

dt 2 


(3.7) 


x n (t) = 

which after taking derivatives leads to 


d n 1 y(t) 
dt 71 - 1 


dx^t) . dy(t) 

= X! = —~~ = X 2 (t) 


dt 
dx-2(t) 
dt 

dx 3 (t) 

dt 


= X 2 = 


= X 3 = 


dt 
d 2 y(t) 
dt 2 

d 3 y(t) 

dt 3 


= x 3 (t.) 

= X 4 (/) 


dx n (t.) _ _ d n y(t) 

— % r. — 


dt 


dy(t) d 2 y(t) 


= -a 0 y(t) - ai ——— - a 2 


dt r 


d n — l 


d n 1 y(t) 


(3.8) 


u(t) 


dt z dt 2 dt 71-1 

= -aox^t) - a 1 x 2 (t) - - a 2 x 3 (t) - -a n _i.x n (/) + u(t) 


The state space form of (3.8) is given by 


X\ 


- 0 

1 

o . 

o - 


" xi(t) 


-o- 
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0 

0 

1 0 ••• 
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(t) 
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0 

0 

. 0 
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Xn— 1 (^) 


0 
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— ai 

-a 2 . 

dn—l - 


. x n (t^ _ 


. 1 . 


u(t) 


(3.9) 
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with the corresponding output equation obtained from (3.7) as 


y(t) = [io 


Xi{t.) 

x 2 (t) 

%n — 1 (1) 

X n (t) 


(3.10) 


The state space form (3.9) and (3.10) is known in the literature as the phase 
variable canonical form. 

In order to extend this technique to the general case defined by (3.5), which 
includes derivatives with respect to the input, we form an auxiliary differential 
equation of (3.5) having the form of (3.6) as 


d n m _ 

dt - + "- 1 


+ a o£(f) = u(t.) 


(3.11) 


for which the change of variables (3.7) is applicable 


.x-i (/) = £(/) 

d£(t) 


X-2 (t) = 
*3(t) = 


dt 

d 2 Z(t) 
dt 2 


(3.12) 


%n{t) 


d^m 

dt 71 - 1 


and then apply the supeiposition principle to (3.5) and (3.11). Since £(/) is the 
response of (3.11), then by the supeiposition property the response of (3.5) is 
given by 


y(t) 


bot(t) + b 1 


d£(t) 

dt 


+ b 2 




(3.13) 


Equations (3.12) produce the state space equations in the form already given 
by (3.9). The output equation can be obtained by eliminating d n £(t.)/dt n from 
(3.13), by using (3.11), that is 


d n m 

dt n 


u(t.) - a n _ix n {t) - - a 1 x 2 (t) - a 0 xi{t) 
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This leads to the output equation 


y(t) = [(60 - «o bn) (h - a t b n ) 


(bn —1 Un—lbn)\ 


'h(i)' 
X-2 ( t) 


+ b n u(t) 


lx n {t)i 

(3.14) 

It is interesting to point out that for b n = 0, which is almost always the case, the 
output equation also has an easy-to-remember form given by 


y(t) = [b 0 bi 


piWl 





(3.15) 


Thus, in summary, for a given dynamic system modeled by differential equation 
(3.5), one is able to write immediately its state space form, given by (3.9) and 
(3.15), just by identifying coefficients a; and bi, i = 0,1,2,...,?? — 1, and using 
them to form the corresponding entries in matrices A and C. 

Example 3.1: Consider a dynamical system represented by the following 
differential equation 


y( 6 ) + 6 y( 5 ) — '2y( 4 ) + y( 2 ) — 5 + 3 y = 7 + u + 4u 


where t/W stands for the /th derivative, i.e. y^ = d l y/df. According to 
(3.9) and (3.14), the state space model of the above system is described by the 
following matrices 


■ 0 

1 

0 

0 

0 

0 ■ 


'O' 

0 

0 

1 

0 

0 

0 


0 

0 

0 

0 

1 

0 

0 

, B = 

0 

0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 


0 

.-3 

5 

-1 

0 

2 

- 6 . 


. 1 . 

c = 

[4 

1 0 

7 

0 

0 ], 

D = 0 
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3.1.2 State Space Variables from Transfer Functions 

In this section, we present two methods, known as direct and parallel program¬ 
ming techniques, which can be used for obtaining state space models from system 
transfer functions. For simplicity, like in the previous subsection, we consider 
only single-input single-output systems. 

The resulting state space models may or may not contain all the modes 
of the original transfer function, where by transfer function modes we mean 
poles of the original transfer function (before zero-pole cancellation, if any, takes 
place). If some zeros and poles in the transfer function are cancelled, then the 
resulting state space model will be of reduced order and the corresponding modes 
will not appeal - in the state space model. This problem of system reducibility 
will be addressed in detail in Chapter 5 after we have introduced the system 
controllability and observability concepts. 

In the following, we first use direct programming techniques to derive the 
state space forms known as the controller canonical form and the observer 
canonical form; then, by the method of parallel programing, the state space 
forms known as modal canonical form and Jordan canonical form are obtained. 

The Direct Programming Technique and Controller Canonical Form 

This technique is convenient in the case when the plant transfer function is 
given in a nonfactorized polynomial form 

Y(s) b n $ n + b n -is n ~ 1 H-F bxs + b 0 

■ 77 — = - 1 - (3.16) 

U(s) s n + cinqs' 11 1 +-h fll-s + a 0 

For this system an auxiliary variable I(s) is introduced such that the transfer 
function is split as 

YM = _ 1 _ 

U( s) s n + a n _is n - 1 +-h ai-s + ciq 

= b n s n + b n _is n 1 + • • • + b\s + bo 

V ( s) 

The block diagram for this decomposition is given in Figure 3.1. 


(3.17a) 

(3.17b) 
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Figure 3.1: Block diagram representation for (3.17) 

Equation (3.17a) has the same structure as (3.6), after the Laplace transforma¬ 
tion is applied, which directly produces the state space system equation identical 
to (3.9). It remains to find matrices for the output equation (3.2). Equation 
(3.17b) can be rewritten as 

T (s) = b n s n \ (s) + b n -\s n H (s) + • • • + b\s\ (s) + boV (s) (3.18) 

indicating that y(t) is just a supeiposition of v(t) and its derivatives. Note that 
(3.17) may be considered as a differential equation in the operator form for zero 
initial conditions, where s = d/dt. In that case, V( 4§, Y( s), and U(s) are simply 
replaced with v(t), //< / ), and u(t), respectively. 

The common procedure for obtaining state space models from transfer 
functions is performed with help of the so-called transfer function simulation 
diagrams. In the case of continuous-time systems, the simulation diagrams 
are elementary analog computers that solve differential equations describing 
systems dynamics. They are composed of integrators, adders, subtracters, and 
multipliers, which are physically realized by using operational amplifiers. In 
addition, function generators are used to generate input signals. The number 
of integrators in a simulation diagram is equal to the order of the differential 
equation under consideration. It is relatively easy to draw (design) a simulation 
diagram. There are many ways to draw a simulation diagram for a given dynamic 
system, and there are also many ways to obtain the state space form from the 
given simulation diagram. 

The simulation diagram for the system (3.17) can be obtained by direct 
programming technique as follows. Take n integrators in cascade and denote their 
inputs, respectively, by v^ n \t), v(t). Use formula (3.18) to 
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construct y(t), i.e. multiply the corresponding inputs r*'/ /) to integrators by the 
corresponding coefficients b; and add them using an adder (see the top half of 
Figure 3.2, where 1/s represents the integrator block). From (3.17a) we have that 

v^ n \t) = u(t) — a n -iv < ' n ~ 1 ^(t) — ••• — — a 0 v(t) 

which can be physically realized by using the corresponding feedback loops in the 
simulation diagram and adding them as shown in the bottom half of Figure 3.2. 



Figure 3.2: Simulation diagram for the direct 
programming technique (controller canonical form) 

A systematic procedure to obtain the state space form from a simulation 
diagram is to choose the outputs of integrators as state variables. Using this 
convention, the state space model for the simulation diagram presented in Figure 
3.2 is obtained in a straightforward way by reading and recording information 
from the simulation diagram, which leads to 
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*(.*.) 


and 


- 0 
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o - 


-o- 
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1 0 ••• 
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0 




0 

x(.f) + 
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0 

O 

O 
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0 

-~a 0 

— a\ 

- a -2 . 

^n—1 - 


.1. 


(3.19) 


y(t) = [(.&„ - «o bn) (h - a t b n ) 


{b n -1 - a n -i6 n )]x(/) + b n u(t) 

(3.20) 


This form of the system model is called the controller canonical form. It 
is identical to the one obtained in the previous section—equations (3.9) and 
(3.14). Controller canonical form plays an important role in control theory since 
it represents the so-called controllable system. System controllability is one of the 
main concepts of modern control theory. It will be studied in detail in Chapter 5. 

It is important to point out that there are many state space forms for a given 
dynamical system, and that all of them are related by linear transformations. 
More about this fact, together with the development of other important state space 
canonical forms, can be found in Kailath (1980; see also similarity transformation 
in Section 3.4). 


Note that the MATLAB function t f 2 s s produces the state space form for 
a given transfer function, in fact, it produces the controller canonical form. 

Example 3.2: The transfer function of the flexible beam from Section 2.6 
is given by 

1.65s 4 - 0.331s 3 - 576s 2 + 90.6s + 19080 
■ ■ “ s 6 + 0.996s 5 + 463s 4 + 97.8s 3 + 12131s 2 + 8.11s 
Using the direct programming technique with formulas (3.19) and (3.20), the 
state space controller canonical form is given by 


'0 

1 

0 

0 

0 

0 ' 


'O' 

0 

0 

1 
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0 

0 


0 

0 

0 

0 

1 

0 

0 

x+ 
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0 

0 

0 

0 

1 

0 

0 

0 

0 

0 

0 

0 

1 


0 

.0 

-8.11 

-12131 

-97.8 

-463 

-0.996. 


.1. 
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and 

y= [19080 90.6 -576 -0.331 1.65 0]x 


o 

Direct Programming Technique and Observer Canonical Form 

In addition to controller canonical form, observer canonical form is related 
to another important concept of modern control theory: system observability. 
Observer canonical form has a very simple structure and represents an observable 
system. The concept of linear system observability will be considered thoroughly 
in Chapter 5. 

Observer canonical form can be derived as follows. Equation (3.16) is written 
in the form 


1' (.s)(.s ri + a n _is n 1 + • • • + a\s + a o) 

= U(s)[b n s n + b n _is n 1 + • • • + b\s + bo) 


(3.21) 


and expressed as 


Y(s) 


-4-(a n _i.s n 1 H-h ais + a 0 )Y(s) 

-\-—U(s)(b n s n + b n _is n 1 + • • • + b\s + bo) 


(3.22) 


leading to 

Y(s) = --a„_iT(.s) - \a n _ 2 Y(s) - -^TT«i^(- s ) - 

S S z s' 1 1 s 


+ b n U(s) + —b n -iU(s) + —^b n - 2 U{ s ) + • • • + b\U{s) + —boU(s) 

(3.23) 

This relationship can be implemented by using a simulation diagram composed 
of n integrators in a cascade, and letting the corresponding signals to pass 
through the specified number of integrators. For example, terms containing 1/s 
should pass through only one integrator, signals a n _ 2 j/(f) and 5„_ 2 "(/) should 
pass through two integrators, and so on. Finally, signals a 0 y(t ) and b 0 u(t ) 
should be integrated //-times, i.e. they must pass through all n integrators. The 
corresponding simulation diagram is given in Figure 3.3. 
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Figure 3.3: Simulation diagram for observer canonical form 


Defining the state variables as the outputs of integrators, and recording 
relationships among state variables and the system output, we get from the above 
figure 

y(t) = x n (t) + b n u(t) (3.24) 

xi(t) = -a 0 y(t) + b 0 u(t) = -a 0 x n (t) + (b 0 - a 0 b n )u(t) 
x 2 (t) = -ary(t) + bruit) + x^t) = x^t) - c^x^t) + ( b n )u(t) 
X 3 (t) = -a 2 y(t.) + b 2 u(t) + x 2 (t) = x 2 (t) - a 2 x n (t) + (b 2 - a 2 b n )u(t) 

X'n(t) — l2/(f) T b n -iU(t) T X n — \(t) 

— $$n — l(f) (In—lXnit) T (bn — 1 —1 b n )u( t ) 

(3.25) 

The matrix form of observer canonical form is easily obtained from (3.24) and 
(3.25) as 
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0 —ao 

0 -cii 

: -a 2 

0 -a n --2 

1 &n — 1 

and 

y(t) = [0 ••• 0 l]x(/) + b n u(t) (3.27) 

Example 3.3: The observer canonical form for the flexible beam from 
Example 3.2 is given by 

'0 0 0 0 0 

1 0 0 0 0 

. _ 0 1 0 0 0 

X “ 0 0 1 0 0 

0 0 0 1 0 

.0 0 0 0 1 

and 

y = [0 0 0 0 0 1 ]x 

o 

Observer canonical form is very useful for computer simulation of linear 
dynamical systems since it allows the effect of the system initial conditions to be 
taken into account. Thus, this form represents an observable system, in the sense 
to be defined in Chapter 5, which means that all state variables have an impact on 
the system output, and vice versa, that from the system output and state equations 
one is able to reconstruct the state variables at any time instant, and of course at 
zero, and thus, determine .i’i( 0 ), .i’ 2 ( 0 )j r w., .r n ( 0 ) in terms of the original initial 
conditions y(0~ ),dy{ 0“ )/dt, ...,d n ~ 1 y(0~ )/dt n ~ 1 . For more details see Section 
5.5 for a subtopic on the observability role in analog computer simulation. 
Parallel Programming Technique 

For this technique we distinguish two cases: distinct real roots and multiple 
real roots of the system transfer function denominator. 


0 19080 

-8.11 90.6 

-12131 , -576 

-97.8 X+ -0.331 “ 
-463 1.65 

-0.996 0 
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Distinct Real Roots 

This state space form is very convenient for applications. Derivation of 
this type of the model stalls with the transfer function in the partial fraction 
expansion form. Let us assume, without loss of generality, that the polynomial 
in the numerator has degree of m < n, then 

W> = _ P m (s) _ 

U(s) (s + Pl )(s + P2 ) • • • (-5 + Pn) 

k i k -2 k n 

= -— + -— + • • • + - — 

S + Pi S + P 2 -5 + Pn 

Here p i. p>- ■■■■p„ are distinct real roots (poles) of the transfer function denom¬ 
inator. 

The simulation diagram of such a form is shown in Figure 3.4. 



Figure 3.4: The simulation diagram for the parallel 
programming technique (modal canonical form) 
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The state space model derived from this simulation diagram is given by 



'-Pi 

0 •• 

. ... o - 


■f 


0 

P'2 ■ ■ 

. ••• 0 


1 

±(t) = 



. 0 

x(f) + 



. 0 

0 •• 

o 

1 

s 3 


.1. 


t(t) 


(3.29) 


y(t) = [&i k- 2 


k n ]x(f) 


This form is known in the literature as the modal canonical form (also known 
as uncoupled form). 

Example 3.4: Find the state space model of a system described by the 
transfer function 


Y{s) _ (s + 5)(-s + 4) 
U(s) (s + 1 )(.s + 2)(.s + 3) 


using both direct and parallel programming techniques. 
The nonfactorized transfer function is 

Y( s) _ s 2 + 9s + 20 
U(s) .s 3 + 6-s 2 + ll.s + 6 


and the state space form obtained by using (3.19) and (3.20) of the direct 
programming technique is 


' 0 

1 

0 ' 


'O' 

0 

0 

1 

x + 

0 

-6 

-11 

— 6_ 


_1_ 


y = [20 9 1 ]x 

Note that the MATLAB function tf2ss produces 


'-6 -11 

-6' 


'1' 

1 0 

0 

X + 

0 

. 0 1 

0 _ 


_0_ 

.. _ N 

n or 




y = [ 1 9 20 ]x 
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which only indicates a permutation in the state space variables, that is 


x = 


0 0 
0 1 
1 0 


1 

0 

0 


x 


Employing the partial fraction expansion (which can be obtained by the 
MATLAB function residue), the transfer function is written as 

Y(s) _ (s + 5)(s + 4) _ 6 6 1 

U(s) (.s + l)(.s + 2)(.s + 3) s + 1 s + 2 + s + 3 

The state space model, directly written using (3.29), is 



'-1 

0 

0 ' 


'1' 

X = 

0 

-2 

0 

x + 

1 


_ 0 

0 

— 3_ 


_ 1 _ 


y = 

[6 

-6 

l]x 



Note that the parallel programming technique presented is valid only for the 
case of real distinct roots. If complex conjugate roots appeal - they should be 
combined in pairs corresponding to the second-order transfer functions, which 
can be independently implemented as demonstrated in the next example. 

Example 3.5: Let a transfer function containing a pair of complex conjugate 
roots be given by 


G(s) = 


S + 1 - j s + 1 + j s + 5 


10 


We first group the complex conjugate poles in a second-order transfer function, 
that is 


G(s) 


85 -[-8 2 3 

.s 2 + 2.S + 2 + s + 5 + s + 10 


Then, distinct real poles are implemented like in the case of parallel programming. 
A second-order transfer function, corresponding to the pair of complex conjugate 
poles, is implemented using direct programming, and added in parallel to the first- 
order transfer functions corresponding to the real poles. The simulation diagram 
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is given in Figure 3.5, where the controller canonical form is used to represent 
a second-order transfer function corresponding to the complex conjugate poles. 
From this simulation diagram we have 


i’l = — 5 .r i + u 
x -2 = —10a; 2 + u 


.X’3 = I4 

■ i ’4 = — 2.1’3 — 2.1-4 + U 
y = 2a; i + 3.12 + 813 + 814 


so that the required state space form is 



'-5 
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0 ' 


' 1 ' 

X = 
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-10 
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X + 
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. 0 

0 

-2 

— 2 _ 


_i_ 


y 

= [2 

3 8 

8 ] 

X 




o 
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Multiple Real Roots 

When the transfer function has multiple real poles, it is not possible to 
represent the system in uncoupled form. Assume that a real pole pi of the transfer 
function has multiplicity r and that the other poles are real and distinct, that is 

Y(s) _ N(s) 

U (S) (S + Pl) T (s + p r+ 1 ) • • • (S + p n ) 

The partial fraction form of the above expression is 

Y ( 5 ) ^11 ^ k\2 k\ r k r ^. 1 k n 

U( s) S + P 1 ( ,S + Pl ) 2 ( s + P 1 ) r S + Pr + 1 s + Pn 

The simulation diagram for such a system is shown in Figure 3.6. 



Figure 3.6: The simulation diagram for the Jordan canonical form 
Taking for the state variables the outputs of integrators, the state space model 
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is obtained as follows 
~~Pi 1 0 

0 -pi 1 0 


0 0 
0 0 


0 ' 
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A = 


0 -pi 

... o 
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-Pi 
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0 

1 

-Pi 

0 

0 
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~Pr +1 


0 

-Pr+2 


0 


B = [0 0 


0 

0 1 1 


• • 1' 


0 

0 

0 

0 

~Pn- 

(3.30) 


C — [^lr klr—1 


^12 ^11 + l ^r + 2 ' ' ' D — 0 


This form of the system model is known as the Jordan canonical form. 
The complete analysis of the Jordan canonical form requires a lot of space and 
time. However, understanding the Jordan form is crucial for correct interpretation 
of system stability, hence in the following chapter, the Jordan form will be 
completely explained. 


Example 3.6: Find the state space model from the transfer function using 
the Jordan canonical form 


G(s) 


s 2 + 6-s + 8 
(s + l) 2 (s + 3) 


This transfer function can be expanded as 


G(s) 


1.25 | 1.5 

•5 + 1 + (S + if 


so that the required state space model is 


0.25 
s T 3 


'-1 

1 

0 ' 


'O' 

0 

-1 

0 

x + 

1 

. 0 

0 

-3. 


_1_ 


y = [1.5 1.25 —0.25]x 


o 
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3.2 Time Response from the State Equation 

The solution of the state space equations (3.1) and (3.2) can be obtained either in 
the time domain by solving the corresponding matrix differential equation directly 
or in the frequency domain by exploiting the power of the Laplace transform. 
Both methods will be presented in this section. 

3.2.1 Time Domain Solution 

For the puipose of solving the state equation (3.1), let us first suppose that the 
system is in the scalar form 


x = ax + bu (3.31) 

with a known initial condition ,r(0) = x 0 . It is very well known from the 
elementary theory of differential equations that the solution of (3.31) is 

t 

x(t) = e at x 0 + j e a ^ t ~ T ' , bu(T)dT (3.32) 

o 

The exponential term e at can be expressed using the Taylor series expansion 
about to = 0 as 


1 \ \ 

e at = 1 + at + ^-a 2 / 2 + + • • • = 'y " Tj( uf )* (3.33) 

~A=0 

Analogously, in the following we prove that the solution of a general nth- 
order matrix state space differential equation (3.1) is given by 

t 

x(/) = e At x(0) + j e A(t - r) Bu(r)r/r (3.34) 

o 

For simplicity, we first consider the homogeneous system without control input, 
that is 


x = Ax, 


x(0) = x, 


(3.35) 
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By analogy with the scalar case, we expect the solution of this differential 
equation to be 

x(f) = e At x(0) (3.36) 


We shall prove that this is indeed a solution if (3.36) satisfies differential 
equation (3.35), where the matrix exponential is defined by using the Taylor series 
expansion as 

11 00 i 

e At = I + At + —A 2 / 2 + —A 3 / 3 + • • • = V —AT* (3.37) 
2! 3! il 

8 = 0 


The proof is simple and is obtained by taking the derivative of the right-hand 
side of (3.37), that is 


de 


At 


d 


dt dt 


= A + — A z t. + -7A 3 / 2 + • • • = A I I + — At. + —A z t. 


1 


= — 1+ A/+ -AN 


2,2 


2 ! 


2 ! 


3! 


2,2 


1! 


2 ! 


= Ae At = e At A 


Now, substitution of (3.36) into differential equation (3.35) yields 
d d 


±= di X= di € X(0) = Ae X(0) = AX(/) 


so that matrix differential equation (3.35) is satisfied, and hence x(f) = < A, x( 0) 
is its solution. 

The matrix e At is known as the state transition matrix because it relates the 
system state at time t to that at time zero, and is denoted by 


$(/) = e At 


(3.38) 


The state transition matrix as a time function depends only on the matrix A. 
Therefore ( I> (/) completely describes the internal behavior of the system, when 
the external influence (control input u( /)) is absent. The system transition matrix 
plays a fundamental role in the theory of linear dynamical systems. In the 
following, we state and verify the main properties of this matrix, which is 
represented in the symbolic form by e At and so far defined only by (3.37). 
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Properties of the State Transition Matrix 

It can be easily verified, by taking the derivative of 


x(.f) = S(.f.)x(0) 

that the state transition matrix satisfies the linear homogeneous state equation 
(3.1) with the initial condition equal to an identity matrix, that is 

= A $(/), $(0) = I (3.39) 

The main properties of the matrix $(/), which follow from (3.37) and (3.38), 
are as follows: 

(a) $(0) = I 

(b) <h _1 (t) = <h(— /) =y <&(/) is nonsingular for every t 

(c) $(/ 2 - to) = $(t 2 ~ h)$(ti - t 0 1 

(d) <!>(/)* = $(■//), for i E iV 

The proofs are straightforward. Property (a) is obtained when / = 0 is 
substituted into the series expansion of <: Al = I + At + • • •. 

Property (b) holds, since 

(e A t )“ 1 e At = I 

which after multiplication from the right by e _At implies 

(e At ) _1 e° = e~ At => $“ 1 (t) = $(-/) 

and (c) follows from 

$(/ 2 - t. 0 ) = e A(t2_t §= eMh-h+u-to) 

= e Mh-t 0 ) e Mn-h) = _ tl )$(/ x _ / Q ) 

Property (d) is proved by using the fact that 

•i>i /! = (e At y = e A(it) = <l>i/7:. 

In addition to properties (a), (b), (c), and (d), we have already established 
one additional property, namely the derivative property, as 
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4-$(/) = A$(/) O ^e At = e At A = Ae At 
at at 

The state transition matrix 4>(/) can be found by using several methods. Two 
of them are given in this chapter—formulas (3.37) and (3.49). The third one, 
very popular in linear algebra, is based on the Cayley-Hamilton theorem and is 
given in Appendix C. 

In the case when the control vector u(f) is present in the system (forced 
response) 

x = Ax + Bu, x(0) = x 0 

we look for the solution of the state space equation in the form 


x(f) = e At i(t) 

Then 

x(f) = Ae At f(f) + e At f(t) = Ax + e At f 
It follows from (3.1) and (3.41) that 

e At f(t) = Bu 


(3.40) 

(3.41) 

(3.42) 


From (3.42) we have 

f(f) = (e At )“ 1 Bu =e“ Al Bu (3.43) 

Integrating this equation, bearing in mind that x(0) = e A '°f(0) = f(0), we get 


f(/) - f(0) = J f" Ar Bu(r)(Ir 
o 

Substitution of the last expression in (3.40) gives the required solution 


(3.44) 


x(f) = e At x(0) + j e A ^ r *Bu(r)r/r 
o 


(3.45) 


or 


x(f) = $(/)x(0) + 


$(/ — t )Bu(r)r/r 


o 


(3.46) 
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When the initial state of the system is known at time t 0 , rather than at time / = 0, 
the solution of the state equation is similarly obtained as 

t 


x(f) = $(/ - / 0 )x(/ 0 ) + 


4>(/ — t )Bu(r)r/r 


to 

t 


= e A(t-to)x( / o)+ j f ,A(i ~ r) Bu(r)flr 


(3.47) 


This can be easily verified by repeating steps (3.40)—(3.45) with x(/ 0 ) = x 0 and 

x(/ 0 ) = e Ato i(t 0 ). 

Example 3.7: For the system given in Example 3.4 find the state transition 
matrix $(/). Evaluate 4>(1) using the MATLAB function expm. Assuming that 
the initial state of the system is zero, find the state response to a unit step. Check 
the solution obtained by using the MATLAB function step. 

At the present time we are able to find the state transition matrix (matrix 
exponential) by using formula (3.37), which deals with an infinite series, and 
hence is not very convenient for calculations. Better ways to find 4>(/) are the 
method based on the Cayley-Hamilton theorem (see Appendix C) and the formula 
based on the Laplace transform, see formula (3.49). However, in this problem, if 
we start with the uncoupled (modal) state space form of the system considered in 
Example 3.4, we can avoid using any of the above methods in order to find the 
state transition matrix. Namely, for diagonal matrices only, it is easy to show that 


-*0 O ' 

0 e~ 2t 0 

0 0 e~ 3t _ 

Using the MATLAB function for evaluating the matrix exponential as 

expm (A* 1), we get 


o 

o 

1 


'0.3679 

0 

0 

0 e“ 2 0 

= 

0 

0.1353 

0 

1 

o 

o 

Cb 

i 

w 

1_ 


0 

0 

0.0498_ 


4>(/) = e L 


-10 0 
0-2 0 
0 0 -3 


4 >( 1 ) 
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The state response to a unit step is computed from (3.46) as 


t 

r 

t 

r 


0 

0 


'i' 

x(/)= / 

4>(f - r)Bu(r)(!r = / 

0 

e-W-T) 

0 


i 

J 

0 

J 

0 

0 

0 

e - 3 h- r ) 


_i_ 


■ e -F-T> ' 


1 — e t 

e -2(t-r) 

dr = 

0.5(1 - e~ 2t ) 

e -3(t~ T ) 


.0.333(1 - e ~ 3t ). 


The step responses of system states, obtained by using MATLAB statements 

[y, x] =step (A, B, C, D) and plot (x) , with 


'-1 

0 

o 


' 1 ' 

0 

-2 

0 

II 

PQ 

1 

1 

o 

0 

1 

CO 

1 


.1. 


are shown in Figure 3.7. 


C = [6 -6 1 ], D = 0 



Figure 3.7: The state responses for Example 3.7 


o 
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3.2.2 Solution Using the Laplace Transform 

The time trajectory of the state vector x(/) can be also found using the Laplace 
transform method. The main properties of the Laplace transform and common 
transform pairs arc given in Appendix A. 

The Laplace transform applied to the state equation (3.1) gives 

sX(s) - x(0“) = AX(s) + BU(s) 
or 

X(s) = (.si - A) _1 x(0") + (.si - A) _1 BU(.s) (3.48) 

Let us assume that x(0) = x(0“). Comparing equations (3.46) and (3.48), it 
is easy to see that the term (si — A) -1 is the Laplace transform of the state 
transition matrix, that is 

4>(.s) = (si - A)” 1 = (|h J Aj adji.sl - A) = C{m} (3-49) 

The time form of the state vector x(f) is obtained by applying the inverse Laplace 
transform to 


X(s) = $(s)x(0) + 4>(s)BU(s) (3.50) 

Note that the second term on the right-hand side corresponds in the time domain 
to the convolution integral, so that we have 

t 


x 


(/) = e At x(0) + / e A(, - r) Bu(r)(Ir 


(3.51) 


Once the state vector x(f) is determined, the output vector y( /) of the system 
is simply obtained by substitution of x(f) into equation (3.2), that is 


y(t) = C4>(/)x(0) + C j 4>(/- t )Bu(r)r/r + Du(t) (3.52) 

o 


or, in the complex domain 


Y(s) = C$(s)x(0) + [C$(s)B + D]U(s) 


(3.53) 
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3.2.3 State Space Model and Transfer Function 

The matrix that establishes a relationship between the output vector Y( s) and the 
input vector U(s), for the zero initial conditions, x(0) = 0, is called the system 
matrix transfer function. From (3.53) it is given by 

G(s) = C(sl - A ) _1 B + D (3.54) 

Note that (3.54) represents the open-loop system matrix transfer function. 

Example 3.8: Find the transfer function for the system given in Example 3.4. 
It is the easiest to use modal (uncoupled) canonical form, which leads to 




's + 1 

0 

0 

' 

-r 

V 


CO 

i 

CO. 

II 

eT 

1] 

0 

s + 2 

0 



1 




0 

0 s + 3j 


_ 1 _ 




r r 

s+1 

0 0 ' 


'l' 




= [6 -6 

1] 

0 

+ 0 


l 






0 

o 


_l_ 







s+3 -1 






6 

6 

_L _ 

1 

4 

+ 5)0 

S + 4) 

” s + 1 

5 + 2 s 

+ 3 (s + 1)( 

s + 2)( 

s + 3) 


If we start with controller canonical form we will get, after some algebra, 


G(s) = [20 9 


s + 6 

-1 0 ' 

-r 

'O' 

11 

s — 1 


0 

6 

0 s . 


_1_ 


= [20 9 


1 


1 

s 3 + 6s 2 + 11s + 6 


s 2 + 6s + 11 
-6 
— 6s 


s + 6 
s(s + 6) 
-ll.s - 6 


1 ' 


'O' 

s 


0 

s 2 


1 


s 2 + 9s + 20 
s 3 + 6s 2 + ll.s + 6 

Note that the MATLAB function ss2tf can be used to solve the above problem. 


o 
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3.3 Discrete-Time Models 

Discrete-time systems are either inherently discrete (e.g. models of bank ac¬ 
counts, national economy growth models, population growth models, digital 
words) or they are obtained as a result of sampling (discretization) of continuous¬ 
time systems. In such kinds of systems, inputs, state space variables, and outputs 
have the discrete form and the system models can be represented in the form 
of transition tables. 

The mathematical model of a discrete-time system can be written in terms 
of a recursive formula by using linear matrix difference equations as 

x[(fc + 1)T] = A d MkT) + B d u(kT) 

(3 55) 

y(kT) = C d x(kT) + B d u(kT) 

Here T represents the sampling interval, which may be omitted for brevity. Even 
more, in the case of inherent discrete systems, there is no need to introduce the 
notion of the sampling interval T so that these systems are described by (3.55) 
with T = 1. 

Similarly to continuous-time systems, discrete state space equations can 
be derived either from difference equations (Subsection 3.3.1) or from discrete 
transfer functions using simulation diagrams (Subsection 3.3.2). In Subsection 

3.3.3 we show how to discretize continuous-time linear systems and obtain 
discrete-time ones. In the rest of the section we parallel most of the results 
obtained in previous sections for continuous-time systems. 


3.3.1 Difference Equations and State Space Form 

An //th-order difference equation is given by 

y(k + n) + a n -iy(k + n - 1) -|- \- a d y(k + 1) + a 0 y(k) 

(3.56) 

= b n u(k + n) + b n _iu(k + n — 1) + • • • + biu(k + 1) + boii(k) 

This equation expresses all values in terms of discrete-time k. 

The corresponding state space equation can be derived by using the same 
techniques as in the continuous-time case. For example, for phase variable 
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canonical form (controller canonical form) in discrete-time, we have 


xi(k + 1) 


‘ 0 1 0 

0 ■ 


X\(k) 


■0‘ 

x -2 (k + 1) 

— 

0 0 1- 

0 


X 2 (k) 

+ 

0 

X n _i(k + 1) 


0 0 0- 

1 


X n — 1(k) 


0 

x n (k+ 1) 


_— (Iq — d i — Cl2 

' &n—l. 


3 


.1. 


y{k) = [(i* 0 - a 0 b n ) ((q - aiM 


(bn—1 dn—lbn)\ 


',l’l ( A?) " 
X 2 (k) 


A-) J 


+ b n u(k) 

(3.57) 

Note that the transformation equations, dual to the continuous-time ones 
(3.11)—(3.13), are given in the discrete-time domain by 


£(k + n) + a n -i£,(k + n - 1) +- \- a^(k + 1) + a 0 £(k) = u(k) (3.58) 

Xl (k) = Z(k) 

X2 (k) = £(k + 1 ) 

x 3 (k) = £(k + 2) ( 3 _ 59 ) 


x n (k) = £(k + n - 1) 

y(k) = b 0 m bi£(k + 1) + b- 2 ^(k + 2) + • • • + b n ^(k + n) (3.60) 

Eliminating ((/,- + n) from (3.60) by using (3.58) and (3.59). the output equation 
given in (3.57) is obtained. 


3.3.2 Discrete Transfer Function and State Space Model 

The derivation of state space equations from discrete transfer functions, based 
on simulation diagrams, is very si mi lar to the continuous-time case. The only 
difference is that in simulation diagrams the integration block 1 /& is replaced by 
the unit delay element z -1 . The state variables are selected as outputs of these 
delay elements. We shall illustrate this method by an example. 
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Example 3.9: Find two state space forms for the transfer function 

y(^) _ s + 1.1 

U(z) ~ U-0.9)(~ + 0.7)(~-0.7) 

We solve this problem by using both direct (a) and parallel (b) programming 
techniques. 

(a) The transfer function can be rewritten as 

y(^) _ s + i.i 

U(z) ~ z 3 - 0.9~ 2 - 0.49. + 0.441 


The simulation diagram for this transfer function is shown in Figure 3.8. 



Figure 3.8: Simulation diagram for direct programming 
in discrete-time domain (controller canonical form) 

The state space model is obtained from this simulation diagram by using the 
outputs of delay elements as state variables. It is given by 



0 

1 

0 ' 


'O' 

x( & + 1) = 

0 

0 

1 

x(/y) + 

0 


0.441 

0.49 

0.9_ 


_ 1 _ 


y(k) = 

[1.1 1 

0]x(&) 



Note that controller canonical form could have been obtained without drawing 
the simulation diagram. We know that this form is identical to phase variable 
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canonical form, which is represented by (3.57). Identifying the corresponding 
coefficients in the original transfer function, the desired state space form is 
obtained directly from (3.57). We have used the above method in order to 
demonstrate at the same time the procedure of drawing simulation diagrams in 
the discrete-time domain. 

(b) Employing the partial fraction expansion (with help of the MATLAB 
function residue), we get 

y(r) _ 6.25 | 0.1786 | -6.4286 

U(z) ~ z — 0.9 + r + 0.7 + r - 0.7 

Since the poles of the transfer function are real and distinct we get the modal 
canonical form as 



'0.9 

0 

0 ' 


'1' 

x(7 + 1) = 

0 

-0.7 

0 

x(7) + 

1 


_ 0 

0 

0.7_ 


_ 1 _ 


y(k) = [6.25 0.1786 -6.4286]x(7) 


o 

3.3.3 Discretization of Continuous-Time Systems 

Real physical dynamic systems are continuous in nature. In this section, we 
show how to obtain discrete-time state space models from continuous-time system 
models. Assume that the plant is linear, continuous, and time invariant with r- 
inputs and p-outputs (see Figure 3.9). Inputs are sampled by using the zero-order 
hold (ZOH) device. This device samples inputs at discrete-time instants kT (see 
Figure 3.10b) and the values obtained for vector u (kT) are held until (k + 1)27 
The corresponding signal is given in Figure 3.10c. 

The state space model of such a plant is 

x(f) = Ax(i) + Bm(i) 

(3.61) 

y(t) = Cx(f) + Dm(i) 

These equations define states and outputs during the sampling interval 
kT < / < (7 + 1)27 Input signals /// i / !. / = l,...,r, are defined by 

nii(t) = m.i(kT) = m(kT ), kT <t < (7+ 1)27 7 = 0,1,2,... (3.62) 
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Figure 3.9: Sampling in a multivariable controlled plant 



Figure 3.10: Transformation of a continuous-time 
input signal by the zero-order hold element 

In the following, we show how to perform discretization of a continuous¬ 
time state space model (3.61) and obtain a discrete-time state space model having 
the form of (3.55) together with the corresponding expressions for matrices 
A d ,B d ,C d , and D d . 
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Consider formula (3.45) with t = T 


x(r) = e AT x(0) + / e A(r “ r) Bu(0)r/r 


o 

T T 

= e AT x(0) + e AT j e“ Ar r/rBu(0) = $(T)x( 0) + j <S>(T - r )drBu(O) 
o o 


which can be written in the form 

x(T) = A d x(0) + B d u(0) 


(3.63) 


(3.64) 


Comparing (3.63) and (3.64) we can find expressions for A,/ and B,/. They are 
given by 

A d = e AT = $(r) 


B,j = e AT / e Ar r/rB = / e A( - T T hlr • B = / e Aa do ■ B 


(3.65) 


0 0 0 

Note that A,/ and B,/ are obtained for the time interval from 0 to T. However, 
it can easily be shown that due to system time invariance the same expressions 
for Ad and B,/ are obtained for any time interval. Namely, steps (3.63)—(3.65) 
can be repeated for succeeding time intervals 2T, 3T,..., (k + 1 )T with initial 
conditions taken, respectively, as x(T), x(2T),.... x( /.' / ). Therefore, for the time 
instant / = (k + 1)T and for t 0 = kT, we have from (3.47) 

(k+l)T 

x[( A- + 1)T] = $[(ifc + 1 )T - kT\x(kT) + I $[(ifc + 1 )T - r]drBu(ifcr) 


(3.66) 


kT 


= A d x(kT) + B d u(kT) 


From the above equation we see that the matrices A,/ and B,/ are given by 
A d = $[( jfc + 1 )T - kT] = $(r) = e AT 

( k + 1 ) T t t (3.67) 

B d = / $[(ifc + 1 )T - r]drB = / <S>(o)doB = I e Aa doB 


kT 


0 


0 
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The last equality is obtained by using change of variables as a = {k + 1)T — r. 
Since (3.65) and (3.67) are identical, we conclude that for a time invariant 
continuous-time linear system, the discretization procedure yields a time invariant 
discrete-time linear - system whose matrices A,/ and B,/ depend only on A, B, 
and the sampling interval T. 

In a similar - manner the output equation (3.61) at / = kT is given by 

y(kT) = Cx(kT) + Du(kT) (3.68) 

Comparing this equation with the general output equation of linear - discrete-time 
systems (3.55), we conclude that 

Cd = C, D d = D (3.69) 


In the literature one can find several methods for discretization of continuous¬ 
time linear - systems. The discretization technique presented in this section is 
known as the integral approximation method. 

In the case of discrete-time linear - systems obtained by sampling continuous¬ 
time linear - systems, the matrix Aj., given by (3.65), can be determined from 
the infinite series 

1 00 i 

A d = e AT = I + AT + —A 2 T 2 + • • • = V -AT (3.70) 

2 ! ?.! 

8 = 0 

It can be also obtained either by using formula (3.49) or the method based on 
the Cayley-Hamilton theorem and setting / = T in <I>(/) = e At . Also, in order 
to evaluate e AT we can use MATLAB function expm(A*T). 

To find B^, the second expression in (3.65) is integrated to give (see 
Appendix C—matrix integrals) 

B d = e AT (-e~ AT A- 1 + A _1 )B = (A d - I)A _1 B (3.71) 


which is valid under the assumption that A is invertible. If A is singular, B^ 
can be determined from 



B,/ = > ' — 


B = 


IT 


r pi J T 1 


A • h! 


A* B 


(3.72) 
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which is obtained by using (3.37) in (3.67) and performing the corresponding 
integration. Note that the above sum converges quite quickly so that only a few 
terms give quite an accurate expression for B d . 

Example 3.10: Find the discrete-time state space model of a continuous¬ 
time system 


x = 


0 1 
-2 -3 


x 


y = [i o ]x 

The sampling period T is equal to 0.1. 

According to (3.65) and (3.69), we have from (3.49) 


a, = $(r) 


- 2e - T -e~ 2T 

e~T - e~ 2T ' 


' 0.9909 

0.0861' 

2e~ 2T - 2e~ T 

■2e- 2T -e- T _ 


-0.1722 

0.7326_ 


B d = (A d -I)A- 1 B 


\Hl + e- 2 T )-e- T ] 


'0.0045' 

e~ T — e~ 2T 


_0.0861_ 


C d = [l 0], D d = 0 

The same result is obtained by using the MATLAB function for discretization of 
a continuous state space model as [ Ad, Bd] = c2d (A, B, T) . 

o 


3.3.4 Solution of the Discrete-Time State Equation 

The objective of this section is to find the solution of the difference state equation 
(3.55) for the given initial state x(0) and the control signal u( k ) at the sampling 
instants T,'2T, ...,kT. For simplicity we assume / = I. 

From the state equation x(k -\- 1) = A d x(k) + B d u(k), for k = 0,1,..., 
N — 1, it follows 

x(l) = A d x(0) + B d u(0) 

x(2) = A d x(l) + B d u(l) = A^x(0) + A d B d u(0) + B d u(l) 

x(3) = A d x('2) + B d u(2) = A^x(0) + A^B d u(0) + A d B d u(l) + B d u(2) 
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N-l 


x(JV) = A d MN - 1) + B d u(N - 1) = A^x(O) + £ A^-^BMi) 


(3.73) 


8 = 0 


Using the notion of the discrete-time state transition matrix defined by 


*d{k) = A d 


(3.74) 


we get 

N -1 

x(JV) = $ d (NM 0) + ^ $ d (N - i - 1 idUum (3.75) 

8 = 0 

Note that the discrete-time state transition matrix relates the state of an input- 
free system at initial time (/,' = 0) to the state of the system at any other time 
k > 0, that is 

x(&) = $d(&)x( 0) = A^x(O) (3.76) 

It is easy to verify that the discrete-time state transition matrix has the 
following properties 

(a) $d(0) = A° = I <= x(0) = $d(0)x(0) 

(b) * d (k 2 - k 0 ) = $ d (k 2 - - k 0 ) = A*’"* 1 A* 1 "*’ = A^~ k ° 

(c) <i> *= (a k d y = Af 

(d) $ d (k + 1) = A d $ d (k), $ d (0) = I 

The last property follows from 

x(k + 1) = A d x(k) =k $ d (k + l)x(0) = A d $ d (k)x(Q) 

It is important to point out that the discrete-time state transition matrix may 
be singular, which follows from the fact that A^ is nonsingular if and only if 
the matrix A d is nonsingular. In the case of inherent discrete-time systems, the 
matrix A d may be singular in general. However, if A d is obtained through the 
discretization procedure of a continuous-time linear system, like in (3.65), then 

(A,)- 1 = (e AT ) _1 = e- Ar 

so that the discrete-time state transition matrix is nonsingular in this case. 
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Remark 1: If the initial value of the state vector is not x(0) but x(& 0 ), 
then the solution (3.75) is 


JV-l 

x (&0 + N) = $ d (N)x.(ko) + ^2 - i - 1 ) B dM k o + *) (3.77) 

8=0 

The output of the system at sampling instant k = N is obtained by substi¬ 
tuting x( /,-) from (3.75) into the output equation, producing 


JV-l 

y (N) = CMNMO) + C d J2^( N ~ i - l.)B d u(*) + D d u (N) (3.78) 

8 = 0 

Note that for / / I, equations (3.75) and (3.78) are modified as 

JV-l 

x(ivr) = $ d (jvr)x(0) + £ ^[(A - * - l)T]B d u(iT) (3.79) 

8 = 0 


JV-l 

y(NT) = CMNTMO) + C d ^ $[( A - * - l)r]B d u(*T) + D d u( AT) 

8 = 0 

(3.80) 

Remark 2: The discrete-time state transition matrix defined by Aj can 
be evaluated efficiently for large values of k by using a method based on the 
Cayley-Hamilton theorem and described in Appendix C. It can be also evaluated 
by using the 2-transform method as given in formula (3.85), see Subsection 3.3.5. 

Example 3.11: For the system given in Example 3.10, use MATLAB to 
find the unit step and impulse responses assuming that the initial condition is 
x(0) = [0 Of. 

The required time responses can be obtained directly by using MAT- 
LAB statements [y,x]=dstep(Ad, Bd, Cd, Dd) (for step response) and 
[y, x] =d±mpulse (Ad, Bd, Cd, Dd) (for impulse response) with the discrete- 
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time model matrices obtained in the last example. The corresponding state and 
output responses are presented in Figure 3.11. 




Figure 3.11: (a) Step responses; (b) impulse responses 


The same problem could have been solved analytically as follows. Since 
the initial condition is zero and u(k) = 1 for k > 0. we get from (3.73) the 
state response as 


JV-l 

x(iV)= £ N = 1,2,... 

8 = 0 

The output response, obtained from (3.78), is given by 


y(N) = Cx(iV), N = 1,2,... 


Flowever, at this point, for large N one is faced with the problem of efficiently 
calculating the powers of matrix A^. This can be facilitated analytically by 
using either the Cayley-Hamilton theorem (see Appendix C) or the Z -transform 
method to be presented in the next subsection. 

By the Cayley-Hamilton method, we have for a 2 x 2 matrix that 

A j = Q'qI + Q'i Ad, k = 2,3,4,... 


— Q'o T 

A 2 = Q'o T Q 1 A 2 


with Q'o and qi satisfying 
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where Ai and A 2 are distinct eigenvalues of A,/. System eigenvalues will be 
considered in Section 3.4. 

o 

3.3.5 Solution of the Discrete State Equation by 
the ^-transform 

Applying the Z -transform (see Appendix B) to the state space equation of a 
discrete-time system 

x(/; + 1) = A d x(k) + B d u(k) (3.81) 

we get 

-X(:) - -:x(0) = A d X(z) + B d U( z) (3.82) 

The complex state vector X(^) can be expressed as 

X(2) = (zl - A d )- 1 ^x(0) + (zl - Adr'BdU(z) (3.83) 

The inverse .r-transform of the last equation gives x( /,•), that is 

x(fc) = Z- 1 [(si - Ad)- 1 ^ ,r(0) + Z- 1 [(si - A d )- 1 B (i U( s)] (3.84) 

Comparing equations (3.75) and (3.84) we conclude that 

$d(k) = 2 _1 [(^I - A^) -1 -:] = A d , k = 1,2,3,... (3.85) 

Let us repeat and emphasize that the discrete state transition matrix <[>,/ (/,•) of 
a general discrete-time invariant linear system can be obtained either by using 
(3.85) or the Cayley-Hamilton method given in Appendix C. 

The inverse transform of the second term on the right-hand side of (3.84) 
is obtained directly by the application of the discrete convolution theorem (see 
Appendix B), leading to 

k — 1 

^{(-I - A (i )- 1 B d U(^)} = ^ $ d (k - 1 - i)B d u(i) (3.86) 

8 = 0 
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Combining (3.84) and (3.86) we get the required solution of the discrete-time 
state space equation as 

k -1 

x(&) = $d(k)M 0) + ^2 ®d(k - i - l)B d u(i) (3.87) 

8 = 0 

The complex form of the output vector Y (s) is obtained if the 2 -transform 
is applied to the output equation 

y{k) = C d x.(k) + D d u(k) 

and X(-:) is substituted from (3.83), leading to 

Y(ff = C d (zl - A d )- 1 ^x(0) + [c d (zl - A^Bd + D d ] U(^) 

From the above expression, for the zero initial condition, i.e. x(0) = 0, the 
discrete matrix transfer function is given by 

G d (z) = C d (zl - Adf'Bd + B d (3.88) 

3.3.6 Response Between Sampling Instants 

An important feature of the state variable method is that it can be modified to 
determine the output between sampling instants. Let to = kT and / = (/,' + A )T, 
where 0 < A < 1. Equation (3.47) gives 

(k+A )T 

x[(jfc + A)T] = C AAT x(kT)+ j e A [(HA)T-r]Bu( r ) ( | r ( 3 . 89 ) 

kT 

Replacing (k + A)T — t by i and assuming that u(r) is constant during 
kT < t < (k + A )T, we get 

AT 

xp + A)T] = e AAT MkT) + 

o 

= AfATMkT) + B d (AT)u(kT) 

where 

A d (AT) = e AAT ( 3 . 91 ) 


e Al3 d/3Bu(kT) 


(3.90) 
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and 


B d (AT) 


AT 


/ 


e Al3 d/3B 


(3.92) 


Therefore, the matrix A d (AT) is obtained by replacing T by AT in A,/. 
Similarly, B^AT) is obtained by replacing T by AT in B^. 


3.3.7 Euler’s Approximation 

Discretization of a continuous-time linear model, as presented in Subsection 3.3.3, 
by the integral approximation method, gives a desired discrete-time linear model. 
However, in the case of high-order systems, computation of the state transition 
matrix is very involved, so that in those cases the matrices A,/ and B,/ are 
calculated approximately by using some simpler methods. The simplest such a 
method, known as Euler’s approximation, is just one of several methods used for 
numerical solution of differential equations. 

The objective of numerical integration is to find a discrete-time counterpart 
to a continuous-time model 


x(f) = A x(f) + Bm(l) 

in the form of a difference equation. The equation obtained, given by a recursive 
formula, is then easily solved by a digital computer. The integration of the 
above equation gives 


x(f) = / [Ax(r) + Bm(r)]r/r 


For simplicity, the main idea of the Euler method is explained for a scalar 
case. Consider the first-order system x = ax + bu. The integration is analogous to 
the problem of finding the area, within the imposed integration limits, between 
the curve defined by f(t) = ax(t) + bu(t) and the time axis. This area is 
approximately equal to the sum of the rectangles in Figure 3.12. 
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Figure 3.12: Euler’s integration method 

If the area is calculated according to Figure 3.12, then from the last expres¬ 
sion for t = (k + 1 )T we get 

kT (k+l)T 


x[(/r + 1)T] = / [Ax(r) + Bm(r)]r/r+ 


[Ax(r) + Bm(r)]r/r 


— OO 


kT 


= x(ifcr) + rAx(l’T) + TBm(kT) 


or 

x[(A: + 1)T] = (I + TA)x(kT) + TBm(kT) (3.93) 

From the last equation, we conclude that for the Euler approximation the state 
and input matrices are given by 

A d = I + T • A, B d = T • B (3.94) 

It can be observed from (3.70), (3.72), and (3.94) that (3.94) produces only 
the first two terms of the series expansion given in (3.70) and only the first term 
of the series expansion given in (3.72). Thus, the Euler approximation is less 
accurate than the integral approximation considered in Subsection 3.3.3, and for 
Euler’s approximation the sampling interval T must be chosen sufficiently small 
in order to get satisfactory results. 

In general, for more accurate computation of the discrete-time model one can 
use any known method for numerical solution of differential equations, e.g. the 
fourth-order Runge-Kutta method or the Adams-Moulton method (Gear, 1971). 
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3.4 The System Characteristic Equation and 
Eigenvalues 

The characteristic equation is very important in the study of both linear time 
invariant continuous and discrete systems. No matter what model type is consid¬ 
ered (ordinary nth-order differential equation, state space or transfer function), 
the characteristic equation always has the same form. 

If we start with a differential nth-order system model in the operator form 

(p n + a n _ip n_1 + • • • + a t p + a 0 )y(t) 

= ( b m p m + + ■■■ + b 1 p + bo)u(t) 

where the operator p is defined as 

: 

p‘ = —-, i = 1,2,...,??.— 1 

and m < n, then the characteristic equation , according to the mathematical 
theory of linear differential equations (Boyce and DiPrima, 1992), is defined by 

s n T (i n —\s n T ■ ■ ■ T r<i'S T f<o = 0 (3.95) 

Note that the operator p is replaced by the complex variable -s playing the role 
of a derivative in the Laplace transform context. 

In the state space variable approach we have seen from (3.54) that 

G(.s) = C(.sl - A) -1 B + D = 1 C[adj(.sl - A)]B + D 

51 

= 1 {C[adj(.sl - A)]B + |,I - A|D} 

5X jTX 

The characteristic equation here is defined by 

| -si — A | = 0 (3.96) 

A third form of the characteristic equation is obtained in the context of the 
transfer function approach. The transfer function of a single-input single-output 
system is 

ru , b m s m + b m -is m 1 + • • • + b\s + bo 

G(s) = -^- 

s n -\- a n -is n + • • • + (i\s -\- (Iq 


(3.97) 
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The characteristic equation in this case is obtained by equating the denominator 
of this expression to zero. Note that for multivariable systems, the characteristic 
polynomial (obtained from the corresponding characteristic equation) appears in 
denominators of all entries of the matrix transfer function. 

No matter what form of the system model is considered, the characteristic 
equation is always the same. This is obvious from (3.95) and (3.97), but is not 
so clear from (3.96). It is left as an exercise to the reader to show that (3.95) 
and (3.96) are identical (Problem 3.30). 

The eigenvalues are defined in linear algebra as scalars, A, satisfying 
(Fraleigh and Beauregard, 1990) 


Av = Av (3.98) 

where the vectors v / 0 are called the eigenvectors. This system of n linear 
algebraic equations (A is fixed) has a solution v 0 if and only if 

| AI — A | = 0 (3.99) 

Obviously, (3.96) and (3.99) have the same form. Since (3.96) = (3.95), it follows 
that the last equation is the characteristic equation, and hence the eigenvalues are 
the zeros of the characteristic equation. For the characteristic equation of order 
n, the number of eigenvalues is equal to n. Thus, the roots of the characteristic 
equation in the state space context are the eigenvalues of the matrix A. These 
roots in the transfer function context are called the system poles, according to the 
mathematical tools for analysis of these systems—the complex variable methods. 
Similarity Transformation 

We have pointed out before that a system modeled by the state space 
technique may have many state space forms. Here, we establish a relationship 
among those state space forms by using a linear - transformation known as the 
similarity transformation. 

For a given system 

x = Ax + Bu, x(0)=x o 
y = Cx + Du 

we can introduce a new state vector x by a linear - coordinate transformation as 
follows 


X = Px 
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where P is some nonsingular n X n matrix. A new state space model is obtained as 

x = Ax + Bu, x(0) = x 0 

(3.100) 

y = Cx + Du 

where 

A = P _1 AP, B = P _1 B, C = CP, D = D, x(0) = P _1 x(0) (3.101) 

This transformation is known in the literature as the similarity transformation. It 
plays an important role in linear control system theory and practice. 

Very important features of this transformation are that under similarity 
transformation both the system eigenvalues and the system transfer function are 
invariant. 

Eigenvalue Invariance 

A new state space model obtained by the similarity transformation does not 
change internal structure of the model, that is, the eigenvalues of the system 
remain the same. This can be shown as follows 

= l.sl- P^APl = IP-^.sI- A)P| 

1 II I (3.102) 

= IP -1 !|.sl — A||P| = |.sl- A| 

Note that in this proof the following properties of the matrix determinant have 
been used 

det(MiM 2 M 3 ) = detMiXdetM 2 xdetM 3 

1 


si - A 


detM -1 = 


detM 


see Appendix C. 


Transfer Function Invariance 

Another important feature of the similarity transformation is that the transfer 
function remains the same for both models, which can be shown as follows 


G(.s) = C .si - A 


-l 


B 


D = CPf.sI - P -1 AP) 1 P _1 B 


D 


= CP[P -1 (sI-A)P] 1 P _1 B + D 


= CPP _1 (.sI - A ) -1 PP -1 B + D 


= C(.sl - A ) -1 B + D = G(.s) 


(3.103) 
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Note that we have used in (3.103) the matrix inversion property (Appendix C) 
(M 1 M 2 M 3)" 1 = M“ 

The above result is quite logical—the system preserves its input-output behavior 
no matter how it is mathematically described. 

Modal Transformation 

One of the most interesting similarity transformations is the one that puts 
matrix A into diagonal form. Assume that P = V = [vi, V 2 ,..., v„], where v; 
are the eigenvectors. We then have 

V -1 AV = A = A = diag {\\, A 2 ,..., A„) (3.104) 

It is easy to show that the elements A;, i = 1,...,??, on the matrix diagonal of 
A are the roots of the characteristic equation |.sl — A| = |-sI — A| = 0, i.e. they 
are the eigenvalues. This can be shown in a straightforward way 

|.si — A| = det{diag(s — Ai, s — A 2 ,..., s — X n )} 

= (s — Ai)(.s — A 2 ) • • • {s — X n ) 

The state transformation (3.104) is known as the modal transformation. 

Note that the pure diagonal state space form defined in (3.104) can be 
obtained only in the following three cases. 

1. The system matrix has distinct eigenvalues, namely Ai^A 2 ^---^A n . 

2. The system matrix is symmetric (see Appendix C). 

3. The system minimal polynomial does not contain multiple eigenvalues. 
For the definition of the minimal polynomial and the corresponding pure 
diagonal Jordan form, see Subsection 4.2.4. 

In the above three cases we say that the system matrix is diagonalizable. 
Remark: Relation (3.104) may be represented in another form, that is 

V -1 A = AV -1 

or 

W r A = AW r 

where 

W T = V 1 =*► W T V = I 
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In this case the left eigenvectors w;. /' = 1.2 . u. can be computed from 

wfA = A;W J =>- A T W; = XiWi 

where W = [wj, w 2 ,..., w n ]. Since |AI —A| = |AI —A r |, then A; is also 
an eigenvalue of A T . 

There are numerous program packages available to compute both the eigen¬ 
values and eigenvectors of a matrix. In MATLAB this is done by using the 
function eig. 

3.4.1 Multiple Eigenvalues 

If the matrix A has multiple eigenvalues, it is possible to transform it into a 
block diagonal form by using the transformation 

J = V _1 AV (3.105) 

where the matrix V is composed of n linearly independent, so-called generalized 
eigenvectors and J is known as the Jordan canonical form. This block diagonal 
form contains simple Jordan blocks on the diagonal. Simple Jordan blocks have 
the given eigenvalue on the main diagonal, ones above the main diagonal with 
all other elements equal to zero. For example, a simple Jordan block of order 
four is given by 

A; 1 0 0 ' 

0 A; 1 0 

0 0 A | 1 

0 0 0 A,-. 

Let the eigenvalue Ai have multiplicity of order 3 in addition to two real 
and distinct eigenvalues, A 2 f A 3 ; then a J matrix of order 5 X 5 may contain 
the following three simple Jordan blocks 

"Ax 1 0 0 O' 

0 Ax 1 0 0 

J = 0 0 Ax 0 0 

0 0 0 A 2 0 

. 0 0 0 0 A 3 . 
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However, other choices are also possible. For example, we may have the 
following distribution of simple Jordan blocks 


-Ai 

1 

0 

0 

0 ' 


‘Ai 

0 

0 

0 

0 ' 

0 

Ai 

0 

0 

0 


0 

Ai 

0 

0 

0 

0 

0 

Ai 

0 

0 

or J = 

0 

0 

Ai 

0 

0 

0 

0 

0 

a 2 

0 


0 

0 

0 

a 2 

0 

o 

0 

0 

0 

A3. 


.0 

0 

0 

0 

A3. 


The study of the Jordan form is quite complex. Much more about the Jordan 
form will be presented in Chapter 4, where we study system stability. 


3.4.2 Modal Decomposition 

Diagonalization of matrix A using transformation x = Vx makes the system 
x = Ax + Bu diagonal, that is 

i = Ak+ (V _1 B)u = Ak+ (W r B)u, x(0) = k 0 

In such a case the homogeneous equation x = Ax, x(0) = x 0 , becomes 

A = Ax, x(0) = V _1 x(0) = V _1 x 0 

or 

•f / — A i-' 1 ' i • i — 1, • • •, n 

This system is represented by n independent differential equations. The modal 
response to the initial condition is 

x(f) = e At x 0 = e At V _1 x 0 = e At W T x 0 


or 

£•;(/) = £q-(0)e Ait = (wfx 0 )e Xlt 

The response x(f) is a combination of the modal components 

x(/) = Vx(/) = V e At V _1 x 0 = Ve At W T x 0 
= (wfx 0 )e Alt V! + (w^x 0 )e A2t v 2 + • • • + (w^x 0 )e A "V 


(3.106) 
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This equation represents the modal decomposition of x( /) and it shows that the 
total response consists of a sum of responses of all individual modes. Note that 
w/xq, are scalars. 

It is customary to call the reciprocals of A; the system time constants and 
denote them by r,-, that is 


Ti = —, i = 1,2,.,.,n 

This has physical meaning since the system dynamics is determined by its time 
constants and these do appeal - in the system response in the form 

The transient response of the system may be influenced differently by 
different modes, depending of the eigenvalues A;. Some modes may decay faster 
than the others. Some modes might be dominant in the system response. These 
cases will be illustrated in Chapter 6. 

Remark: A similarity transformation A = V -1 AV can be used for the 
state transition matrix calculation. Recall 

x(/) = e At x(0), x(f) = Vx(f), x(0) = Vx(0) 


and 

x(/) = V _1 e A ‘Vx(0) = $(/)x(0) 


Hence, 

$(/) = e At = V _1 e At V = W T e At V 
or, in the complex domain 


$(s) = V -1 (.sI - A) _1 V 
= V~ 1 diag{s — Ai, s — X- 2 ,..., s — A n } _1 V 
1 1 1 


= V 1 diag 


s — Ai .s — A '2 


s — \ r 


V 


Remark: The presented theory about the system characteristic equation, 
eigenvalues, eigenvectors, similarity and modal transformations can be applied 
directly to discrete-time linear systems with A,/ replacing A. 
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3.5 State Space MATLAB Laboratory Experiments 

In this section we present three MATLAB laboratory experiments on the state 
space method in control systems. These experiments can be used either as sup¬ 
plements for lectures or independently in the corresponding control system lab¬ 
oratory. Most of the required MATLAB functions have been already introduced 
in the examples done in this chapter. Students should also consult Appendix 
D, where a shortened MATLAB manual is given. It is advisable that before 
using any MATLAB function, the students check all its options by typing help 
function name. 


3.5.1 Experiment 1—The Inverted Pendulum 

Part 1. The linearized equations of the inverted pendulum, obtained by 
assuming that the pendulum mass is concentrated at its center of gravity (Kwak- 
ernaak and Sivan, 1972; Kamen, 1990) are given by 


(./ + m,L 2 )0(t) — mgL0(t) + mLd(t) = 0 
(M + m)d(t) + mLQ(t) = u(t) 


(3.107) 


where 9(t) is the angle of the pendulum from the vertical position, d(t) is the 
position of the cart. u(t) is the force applied to the cart. M is the mass of the 
cart, m is the mass of the pendulum, g is the gravitational constant, and J is the 
moment of inertia about the center of mass. Assuming that normalized values are 
given by J = L, /. = % g = 9.81, M = 1, m = 0.1, derive the state space form 


x(f) = A x(f) + B u(t) 


where 

x(f) = [0(t) 9(t) d(t) d(t)] T 

and A 4x4 and B 4xl are the corresponding matrices. 

Part 2. Using MATLAB determine the following: 

(a) The eigenvalues, eigenvectors, and characteristic polynomial of matrix A. 

(b) The state transition matrix at the time instant / = I. 

(c) The unit impulse response (take 9(t) and d(t) as the output variables) for 
0 < / < 1 with the step size At = 0.1 and draw the system response using 
the MATLAB function plot. 
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(d) The unit step response for 0 < / < 1 and At = 0.1. Draw the system 
response. 

(e) The unit ramp response for 0 < / < I and At = 0.1 and draw the system 
response. Compare the response diagrams obtained in (c), (d), and (e). 

(f) The state response resulting from the initial state x(0) = [—1 1 1 I ] T 

and the input u(/) = sin(f) for 0 < / < 5 and At = 0.1. 

(g) The inverse of the state transition matrix (e At ) 1 for t = 5. 

(h) The state x(f) at time / = 5 assuming that x(10) = [10 0 5 2] r and 

u(t) = 0 by using the result from (g). 

(i) Find the system transfer function. 

Part 3. Discretize the continuous-time system given in (3.107) with 
T = 0.02 and find the discrete space model 

x(/r + 1) = A d -x(k) + B d u(k) 


Assuming that the output equation of the discrete system is given by 


y(A) 


1 o o 
o o 1 


x(&) = C d x.(k) 


find the system (output) response for 0 < /.■ < 50 due to initial conditions 
x 0 = [—1 1 —1 I ] T and unit step input (note that u(k) should be generated 

as a column vector of 50 elements equal to 1). 

Part 4. Consider the continuous-time system given by 


d 2 y(t) 

dt 2 


+ 0.1 


dy(t) 

dt 


u(t) 


(3.108) 


(a) Discretize this system with T = 1 by using the Euler approximation. 

(b) Find the response of the obtained discrete system for k = 1, 2,3,...,20, when 
u(t) = sin(0.l7r/) and j/(0) = y(0) = 0. 

(c) Find discrete transfer function, characteristic equation, eigenvalues, and 
eigenvectors. 

Part 5. Discretize the state space form of (3.108) obtained by using 
MATLAB function c2d with / = I. Find the discrete system response for 
the initial condition and the input function defined in Part 4b. Compare the 
results obtained in Parts 4 and 5. Comment on the results obtained. 
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3.5.2 Experiment 2—Response of Continuous Systems 

Part 1. Consider a continuous-time linear system represented by its transfer 
function 

G{S) = .s 2 + 5 s + 6 

(a) Find the impulse response by using the MATLAB function impulse. In this 
case you have to use [y, x] =±mpulse (num, den) , where num and den 
are row vectors that contain the polynomial coefficients in descending powers 
of s. Plot both state space variable and output responses (use function plot). 

(b) Find the step response by using the function step and plot both the state 
response and the output response. 

(c) Find the zero-state response due to an input given by f(t) = e 3t , / > 0. 
Note that you have to use the function lsim and specify input at every 
time instant of interest. That can be obtained by t=0 : 0.1: 5 (defines t 
at ().(). I. 0.2,.,... I . 9 , 5 ) and by f = exp(-3*t). Check that the results 
obtained in (c) agree with analytical results at t = 1. 

(d) Obtain the state space form for this system by using the function tf2ss. 
Repeat parts (a), (b), and (c) for the corresponding state space representation. 
Use the following MATLAB instructions 

[y,x]=impulse(A,B,C,D,1); 

[y,x]=step(A,B,C,D,1) ; 

[y, x] = lsim (A, B, C, D, f, t) ; 

respectively, with / and t defined in (c). Compare the results obtained. 
Part 2. Consider the continuous-time linear - system represented by 

f(t) = e~ 4 \ t > 0, y{ 0") = 2, y( 0) = 1 

(a) Find the complete system response by using the MATLAB function lsim. 
Compare the simulation results obtained with analytical results. Hint: Use 

[y, x] =lsim (A, B, C, D, f, t, X0) ; 

with t = 0:0.1:5. Note that the initial condition for the state vector, X0, has 
to be found. This can be obtained by playing algebra with the state and 
output equations and setting / = 0. 
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(b) Find the zeros and poles of this system by using the function tf2zp. 

(c) Find the system response due to initial conditions specified in Paid 2a and 
the impulse delta function as an input. Since you are not able to specify the 
system input in time (the delta function has no time structure), you cannot use 
the lsim function. Instead use the initial function (zero-input response). 
The required response is obtained analytically as follows 

XI / ! = < A I XI 0 ! + B) 

where A and B stand for the system and input matrices in the state space. 
Thus, the new initial condition is given by x(0) + B. 

(d) Justify the answer obtained in (c). Solve the same problem analytically by 
using the Laplace transform. Plot results from (c) and compare with results 
obtained in (d). Can you draw any conclusion for this “nonstandard” problem 
from the point of view of the system initial conditions at / = 0 + . (The 
standard problem requires that for the impulse response all initial conditions 
are set to zero.) 

Part 3. Given the following dynamical system represented in the state space 
form by (Gajic and Shen, 1993) 


' -0.01357 

-32.2 

-46.3 

0 ' 


' -0.433 ' 

0.00012 

0 

1.214 

0 

, B = 

0.1394 

-0.0001212 

0 

-1.214 

1 

-0.1394 

0.00057 

0 

-9.1 

-0.6696. 


-0.1577. 


C = 


'o o o r 
.1 o o oj’ 


D = 0 2xl 


This is a real mathematical model of an F-8 aircraft (Teneketzis and Sandell, 
1977). Using MATLAB, determine the following quantities. 


(a) The eigenvalues, eigenvectors, and characteristic polynomial. Take 
p=poly (A) and verify that roots (p) produces also the eigenvalues 
of matrix A. 

(b) The state transition matrix at the time instant / = 1. Use the expm function. 

(c) The unit impulse response and plot output variables. Hint: Use 

impulse(A, B, C, D) ; 
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(d) The unit step response and plot the corresponding output variables. 

(e) Let the initial system state be x(0) = [—1 1 0.5 1 ] T . Find the response 

due to an input given by /(/) = sin(f), 0 < / < 1000. Hint: Take 
t=0:10:1000 and find the corresponding values for f(t) by using the 
function sin in the form / = sin(f). Then use the lsim function. 

(f) Find the system transfer functions. Note that you have one input and two 
outputs which implies two transfer functions. Hint: Use the function ss2tf. 

(g) Find the inverse of the state transition matrix (e At ) = e _At at / = 2. 

Part 4. Consider a linear continuous-time dynamical system represented by 

its transfer function 

r , (s + l)(-s + 3 )(s + 5 )(.s + 7) 

' ' .s(.s + 2)(.s + 4)(.s + 6)(.s + 8)(.s + 10) 

(a) Input the system zeros and poles as column vectors. Note that in this case 
the static gain k = 1. Use the function zp2ss (z, p, k) in order to get the 
state space matrices. 

(b) Find the eigenvalues and eigenvectors of matrix A. 

(c) Verify that the transformation x = Px, where P is the matrix whose columns 
are the eigenvectors of matrix A. produces in the new coordinates the 
diagonal system matrix A = P _1 AP with diagonal elements equal to the 
eigenvalues of matrix A. 

(d) Find the remaining state space matrices in the new coordinates. Find the 
transfer function in the new coordinates and compare it with the original 
one. 

(e) Compare the unit step responses of the original and transformed systems. 

3.5.3 Experiment 3—Response of Discrete Systems 

Part 1. Consider a discrete-time linear system represented by its transfer 

function 

(a) Find the impulse response by using the MATLAB function dimpulse. In 
this case you have to use [y, x] =d±mpulse (num, den), where num and 
den are row vectors which contain the polynomial coefficients in descending 
powers of -s. Plot both state and output responses (use function plot). 
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(b) Find the step response by using the function dstep and plot both the state 
and output responses. 

(c) Find the system (output) response due to a unit step function, f(k) = h(k), 
and initial conditions specified by y( — 1) = 0, y(—'2) = 1. Note that you 
have to use the function dlsim and to specify input at every time instant of 
interest. That can be obtained by k=0 :1: 20 (defines k at 0,1,2,..., 19,20) 
and by f (k) =1. Check analytically that the results obtained in (c) agree 
with the analytical results for k = 10. 

(d) Obtain the state space form for this system by using the function tf2ss. 
Repeat parts (a), (b), and (c). Use the following MATLAB statements 

[y,x]=d±mpulse(A,B,C,D,1); 

[y,x]=dstep(A,B,C,D,1); 

[y, x] =dls±m (A,B,C,D,f,k); 

respectively, with / and k defined in (c). Compare the results obtained. 
Part 2. Consider the discrete-time linear system represented by 

y(k + 2) + ^y(k + 1) + ^(/r) = f(k + 1) 
f(k)=(0.8) k u(k ) i:; y(- 1) = 2, y(- 2) = 3 

(a) Find the system response by using the MATLAB function dlsim. Hint: Use 

[y,x]=dlsim(A,B,C,D,f,X0) ; with k=0 :1:10. 

Note that the initial condition has to be found. This can be obtained by 
playing algebra with the state space and output equations. Compare the 
simulation results obtained with analytical results. 

(b) Find the zeros and poles of this system by using the function tf2zp. 

(c) Find the system response due to initial conditions specified in Paid 2a and 
with the impulse delta function as an input. Use the dlsim function. 

(d) Solve the same problem analytically by using the Z -transform. Plot results 
from (c) and compare with results obtained in (d). 

Part 3. Given a dynamical system represented in the continuous-time state 
space form in Section 3.5.2, Experiment 2, Paid 3. 

(a) Discretize the continuous-time system by using the MATLAB function c2d. 
Assume that the sampling period is T = 1. 
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(b) Find the eigenvalues, eigenvectors, and characteristic polynomial of the 
obtained discrete-time system. 

(c) Find the state transition matrix at time instant k = 5. 

(d) Find the unit impulse response and plot output variables. Flint: Use 

dimpulse(A,B,C,D); 

(e) Find the unit step response and plot the corresponding output variables. 

(f) Assume that the initial system state is x(0) = [— 1 0 1 — 0.5] T . Find 
the response due to an input given by f(t) = sink, 0 < k < 1000. Flint: 
Take k=0 :10 :1000 and find the corresponding values for f(k) by using 
the function sin in the form / = sin(fc). Then use the dlsim function. 
Compare the obtained discrete-time results with the continuous-time results 
for the same system studied in Section 3.5.2, Experiment 2. 

(g) Find the system transfer functions. Note that you have one input and two 
outputs which implies two transfer functions. The matrices C and D are not 
changed due to discretization procedure. Hint: Use the function ss2tf. 
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3.7 Problems 

3.1 An antenna control problem (Dressier and Tabak, 1971) is represented by 
the open-loop transfer function 

= _ A (- s + 1 ) _ 

s 2 (s + 6 )(s + 11.5)(.s' 2 + 8 s + 256) 

Find state space matrices for the following forms: 


(a) Controller canonical form. 
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(b) Observer canonical form. 

3.2 A robotic manipulator called the acrobot has the following linearized model 
(Spong. 1995) 


' 0 

0 

1 

O' 


' 0 ' 

0 

0 

0 

1 

, B = 

0 

12.49 

-12.54 

0 

0 

-2.98 

.-14.49 

29.36 

0 

0. 


. 5.98 . 


Assume that the output matrices are given by 

C = [1 0 1 0], D = 0 


Use MATLAB in order to find the following quantities: 

(a) Eigenvalues and characteristic polynomial. 

(b) Modal canonical form. 

(c) Open-loop transfer function. 

(d) Controller and observer canonical forms. 

3.3 Consider the harmonic oscillator in the state space form 



0 1 .i’i 

-1 0 X-2 



y = [o 



•i'i(O) 

* 2 ( 0 .). 



(a) Find the state transition matrix. 

(b) Find the system response due to a unit step input. 

(c) Verify the answer obtained by using the MATFAB functions ss2zp 
and lsim. 

3.4 Given the matrix 



Find e At by the Cayley-Hamilton and the Faplace transform methods. 

3.5 For the system 


y(t) + 2 y(t) + y(t) = Qii(t) + u(t) 
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(a) Draw the simulation block diagram. Use any method. 

(b) Find the system impulse response by the MATLAB function im¬ 
pulse. 


3.6 


Given a discrete system x(k + 1) = A xi /.•),• where 

'1 2 ' 


A = 


0 3 


Find its response due to the initial condition given by x(0) = [1 
3.7 Given a linear time invariant continuous system 

ro' 


ir- 


x(/) = Ax(t) 


1 


u(t) 


with 


x 


a) = 


Xl(l) 

1 * 2 ( 1 .)] 


$ " i/! = { 


1 / > 2 
0 0 < / < 2 


Assuming that the state transition matrix has a known (given) form as 

ft 11 '* - ' 0 ! ° ' 

\_4 > 2l{t ~ to) </>22(t — to) . 

find the system response for any / > 0. 

3.8 Find the impulse response of the system 

y + '2y + 10 y = u — 3ii + 5 u 


Find the system transfer function and the state space form. 

3.9 Find the system response for / > I due to its initial condition at / = I 

+ I//!/! + 3 J y(T)dT = 0, 2/(1) = 2 


3.10 Find the response of the discrete system 

y(k + -2)+ y(k)= (-1)*, 2 /( 0 )= 1, 2 /( 1 ) = 0 

Verify the answer by the MATLAB function dlsim. 

3.11 A continuous-time system is represented by 

y + 4y + 3y = u 

(a) Find the transfer function and the impulse response. 

(b) Compute the response y(t) for y(0 ~) = —2, y(0~) = 1 and u(t) 
equal to a unit step function. 
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3.12 Given a linear continuous-time system with 



' 0 1 


' 0 


'-1' 

A = 

-2 -3 

,B = 

-1 

, x(0) = 

0 


(a) Find the state transition matrix. 

(b) Find the system transfer function. 

(c) Find the system response due to a unit step input. 

(d) Verify the answers obtained by using MATLAB. 

3.13 A discrete system is given by 

y(k + 1) — 0.5 y(k) = 2 u(k + 1) + u(k) 


Compute the impulse response. Verify the result by the MATLAB function 

dimpulse. 

3.14 Discretize the following system by using the Euler approximation 

y + 2y + 3sin (y(t)) = u(t), y( 0) = 1, y( 0) = 2 


3.15 Given a time invariant linear system with the impulse response equal to 
e~*. Find the response of this system due to an input '26(t — 1) + 3 h(t — 2). 
where 6(t) is the impulse delta function and h(t) is a unit step function. 
What MATLAB functions can be used to solve this problem? 

3.16 A lineal - discrete system is represented by 


A = 


0 

-2 



B = 



, C = [0 1], D = 0 


(a) Find its state transition matrix. 

(b) Find the transfer function. 

(c) Find the system response due to u(k) = k assuming that .i’i(0) = 1 
and .r 2 (0) = 3. 


3.17 Find the response of the system given by 

y+2y + y = ii + u, u(t) = 2e _t , r/(0) = 1, r/(0) = 1 
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3.18 Given a second-order linear system at rest (initial conditions are zero) 

y+ ' 2 &nV + uly = ulu(t) 

Find its unit step response for ( < I. 

3.19 Find the response of the discrete system 

y(k + 2) - 6 y(k + 1) + 8 y(k) = 3 k + 2, 2 /( 0 ) = 1, y( 1) = 1 

3.20 Find the response of the continuous system 

y + 3y-10y = 2u + bu, y(0) = 1, y(0) = -1, u(t) = t 

3.21 Discretize the system y = u by using both the Euler and the integral 
approximations. Compare the discrete systems obtained. 

3.22 Given a linear continuous system 

x(f) = A x(f) 


with 

Find the similarity transformation such that this system has the diagonal 
form in the new coordinates. 

3.23 Find the state transition matrix of a continuous system with 



Use the Taylor series expansion method. 

3.24 Find the response of a discrete system represented by 

y(k + 2) + 2 y(k) + 1 = (- l) fc , 2/(0) = J/G) = 1 




156 


STATE SPACE APPROACH 


3.26 Consider a fifth-order industrial reactor (Arkun and Ramakrishnan, 1983; 
Petkovski et ah, 1991) represented by 


"—16.11 

-0.39 

27.2 

0 

0 

0.01 

-16.99 

0 

0 

12.47 

15.11 

0 

-53.6 

-16.57 

71.78 

-53.36 

0 

0 

-107.2 

232.11 

2.27 

60.1 

0 

2.273 

-102.99 

'11.12 

-2.61 

-21.91 

-53.5 

69.ll T 

-12.6 

3.36 

0 

0 

0 


0 0 0 0 1 
“ [o 1 1 0 0_ 

Using MATLAB, find the following: 

(a) The system transfer function. 

(b) The impulse response. 

(c) The response due to inputs u\{t) = e~ t + sin(t) and u- 2 (t) = 0. 


3.27 Discretize the system given in Problem 3.26 by the MATLAB function c2d 
with T = 0.1 and repeat the steps (a), (b), and (c) from Problem 3.26. 

3.28 The model of a synchronous machine connected to an infinite bus (Koko- 
tovic et ah, 1980; Grodt and Gajic, 1988) has the system matrix 


'-0.58 

0 

0 

-0.27 

0 

0.2 

0 

0 

-1 

0 

0 

0 

1 

0 

0 

0 

-5 

2.1 

0 

0 

0 

0 

0 

0 

0 

337 

0 

0 

-0.14 

0 

0.14 

-0.2 

-0.28 

0 

0 

0 

0 

0 

0 

0 

0.08 

2 

.-17.2 

66.7 

-11.6 

40.9 

0 

-66.7 

-16.7 


(a) Find the eigenvalues, eigenvectors, and similarity transformation that 
puts this system into a diagonal form. 

(b) Discretize this system with 7 = I. 

(c) Find the response of the discrete system obtained in part (b) due 

to the initial condition x(0) = [1 1 1 1 1 1 1] T and draw 

the corresponding response for the time interval 0 < k < 10. Use 
MATLAB. 
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3.29 A linearized mathematical model of an aircraft considered in Litkouhi 
(1983) and Gajic and Shen (1991) has the form 


-0.015 

-0.0805 

-0.0011666 

0 

0 

0 

0 

0.03333 

-2.28 

0 

-0.84 

1 

0.6 

0 

-4.8 

-0.49 


-0.0000916 0.0007416 
0 0 

- 0.11 0 

-8.7 0 



Obtain the following (using MATLAB): 

(a) Discretize this model with 7 = I. 

(b) Find its response due to a unit ramp input. 

(c) Find the system transfer function and the system poles. 

3.30 Show by induction that the characteristic equation (3.96) of a system in the 
phase variable canonical form is indeed given by (3.95). 

3.31 Write a MATLAB program to obtain the system’s modal form for Example 
3.4. Using that program, check the corresponding results in Examples 3.4 
and 3.9. 



